clear
version 14
set more off
cap log close

mata: mata clear
mat drop _all

gl OUT 	"${LOCDRIVE}\Results"
gl TEMP	"${LOCDRIVE}\Tempdata"
gl FIG "${LOCDRIVE}\Figures"

loc r1_hh hhid
loc r1_weight hh_weight
loc r1_indiv sbmemno
loc r2_hh y2_hhid
loc r2_weight y2_weight
loc r2_indiv indidy2
loc r3_hh y3_hhid
loc r3_weight y3_weight
loc r3_indiv indidy3

loc ea ea
loc strata strataid
loc admin0 macroregion
loc admin1 region 		
loc admin2 district  	
loc admin3 ward

loc numeraire other

loc cre 1 /*include correlated random effects*/
loc crelev hh_c /*level of correlated random effects*/ 

loc servicelab "Services"
loc food1lab "SF"
loc food2lab "PNFFV"
loc food3lab "ASF"
loc food4lab "OF"
loc foodlab "Food"
loc goodlab "Goods"
loc otherlab "Other"

loc bootreps=10 /*Note the manuscript uses 300 bootreps, set low for testing purposes*/
loc maxiter=300
loc tot_didnotconverge = 0

forvalues r = 1/3 {
	loc xc_dfl_r`r' = ${XC_DFL_R`r'}
	loc dfl_r`r' = ${DFL_R`r'}
	}

log using "${TEMP}\Analysis_SE_Boot_SL_food${p1_endog}good${p2_endog}service${p3_endog}_np${npowers}", replace

display "$endogvarlist"
display "$exogvarlist"
display "$eqlist"

** Can change these if needed - they are set in the previous estimation file
forvalues i=1/3 {
	display "`i' endog: $p`i'_endog" 
	}


** First draw samples from urban and rural strata

use "${TEMP}/Panel_all_estready_3sls_subgroups.dta", clear
merge 1:1 `r1_hh' round using "${TEMP}\Tanzania_Panel_estready_all.dta", nogen keep(1 3) keepusing(admin3_c)

tempfile sample_0
save `sample_0', replace

g count = 1
keep if round==1
collapse (sum) count, by(admin3_c urb rur) 
count if urb==1
loc bootsize_urb = r(N)-1
count if rur==1
loc bootsize_rur = r(N)-1

foreach pop in urb rur {
    use `sample_0', clear
	keep if round==1 & `pop'==1  
	
	keep `r1_hh' admin3_c urb rur 

	tempfile sample_list_`pop'
	save `sample_list_`pop'', replace
	}

	
set seed 210627
forvalues b=1/`bootreps' {
    display as error "Boot rep is: `b'/`bootreps'"
	
	*time stamp
	local today =  "`c(current_date)' `c(current_time)'"
	display as error "Time stamp: `today'"

	foreach pop in urb rur {
		use `sample_list_`pop'', clear
		bsample `bootsize_`pop'', cluster(admin3_c) 
		tempfile samp_`pop'
		save `samp_`pop'', replace
		}
			
	use `samp_rur', clear
	append using `samp_urb'
	
	g obs = _n
	expand 3
	g count=1
	bysort obs: g round = sum(count)
	drop count obs
	
	merge m:1 `r1_hh' round using `sample_0', nogen keep(1 3) 
	
	replace y=y_stone
	g y_backup=y_stone
	g y_old=y_stone
	g y_change=0
	scalar crit_test=1
	scalar iter=0
	
	while crit_test>$conv_crit & iter<`maxiter' {
	scalar iter=iter+1
	quietly reg3 $eqlist [aweight=wt_cur], constr($conlist) endog($endogvarlist) exog($exogvarlist)
	if (iter>1) {		
		matrix params_old=params
		}
	matrix params=e(b)
	quietly replace pAp=0
	quietly replace pBp=0
	quietly replace y_old=y
	quietly replace y_backup=y

	*predict with y=1
	*generate rhs vars,interactions with y=1
	forvalues j=1(1)$npowers {
		quietly replace y`j'=1
		} 
	forvalues j=1(1)$neq {
		quietly replace ynp`j'=np`j'
		tempvar crenpy 											
		bysort `crelev': egen `crenpy' = mean(np`j')
		quietly replace crepy`j'=`crenpy'
		} 
	forvalues j=1(1)$ndem {
		quietly replace yz`j'=z`j'
		} 
		
	*generate predicted values
	forvalues j=1(1)$neq {
		quietly predict s`j'hat_y1, equation(s`j')		
		}
		
	*set all p, pz, py to zero
	foreach yvar in $nplist $ynplist $npzlist $crepylist $creplist {  
		quietly replace `yvar'=0
		} 
	forvalues j=1(1)$neq {
		quietly predict s`j'hat_y1_p0, equation(s`j')		
		}
	
	*refresh p,pz
	forvalues j=1(1)$neq {
		quietly replace np`j'=np`j'_backup
		tempvar crenp											
		bysort `crelev': egen `crenp' = mean(np`j'_backup)
		quietly replace crep`j' = `crenp'
		forvalues k=1(1)$ndem {
			quietly replace np`j'z`k'=np`j'_backup*z`k'
			}
		}
	
	*generate rhs vars,interactions with y=0
	foreach yvar in $ylist $ynplist $yzlist $crepylist { 
		quietly replace `yvar'=0
		} 
		
	*generate predicted values
	forvalues j=1(1)$neq {
		quietly predict s`j'hat_y0, equation(s`j')		
		}
		
	*set all p, pz, py to zero
	foreach yvar in $nplist $ynplist $npzlist $crepylist $creplist { 
		quietly replace `yvar'=0
		} 
	forvalues j=1(1)$neq {
		quietly predict s`j'hat_y0_p0, equation(s`j')		
		}

	*refresh p only
	forvalues j=1(1)$neq {
		quietly replace np`j'=np`j'_backup
		}
	
	*fill in pAp and pBp
	forvalues j=1(1)$neq {
		quietly replace Ap`j'=s`j'hat_y0-s`j'hat_y0_p0
		quietly replace pAp=pAp+np`j'*Ap`j'	
		quietly replace Bp`j'=(s`j'hat_y1-s`j'hat_y1_p0)-(s`j'hat_y0-s`j'hat_y0_p0)
		quietly replace pBp=pBp+np`j'*Bp`j'	
		quietly drop s`j'hat_y0 s`j'hat_y0_p0 s`j'hat_y1 s`j'hat_y1_p0
		}

	*round pAp and pBp to the nearest millionth, for easier checking
	quietly replace pAp=int(1000000*pAp+0.5)/1000000
	quietly replace pBp=int(1000000*pBp+0.5)/1000000

	*recalculate y,yz,py,pz
	quietly replace y=(y_stone+0.5*pAp)/(1-0.5*pBp)
	forvalues j=1(1)$npowers {
		quietly replace y`j'=y^`j'
		}
	forvalues j=1(1)$ndem {
		quietly replace yz`j'=y*z`j'
		} 

	*refresh py,pz
	forvalues j=1(1)$neq {
		quietly replace ynp`j'=y*np`j'_backup
		tempvar crenpy 										
		bysort `crelev': egen `crenpy' = mean(ynp`j')
		quietly replace crepy`j'=`crenpy'
		tempvar crenp 										
		bysort `crelev': egen `crenp' = mean(np`j'_backup)
		quietly replace crep`j'=`crenp'
		forvalues k=1(1)$ndem {
			quietly replace np`j'z`k'=np`j'_backup*z`k'		 
			}
		}

	if (iter>1 & conv_param==1) {		
		matrix params_change=(params-params_old)
		matrix crit_test_mat=(params_change*(params_change'))
		svmat crit_test_mat, names(temp)
		scalar crit_test=temp
		drop temp
		}
	quietly replace y_change=abs(y-y_old)
	quietly summ y_change
	if(conv_y==1) {
		scalar crit_test=r(max)
		}
	display "iteration " iter 
	scalar list crit_test 
	summ y_change y_stone y y_old pAp pBp
	}

	loc didnotconverge = (iter>=`maxiter')
	
	if `didnotconverge'==0 { /*proceed if it converged on the prior step*/

		*now, create the instrument, and its interactions yp and yz
		quietly replace y_inst=(y_tilda+0.5*pAp)/(1-0.5*pBp)
		forvalues j=1(1)$npowers {
			quietly replace y`j'_inst=y_inst^`j'
			}
		forvalues j=1(1)$nprice { 
			if ${p`j'_endog}==1 {
				replace ynp`j'_inst1=y_inst*p`j'_inst1
				replace ynp`j'_inst2=y_inst*p`j'_inst2	
				replace ynp`j'_inst3=y_inst*p`j'_inst3	
				}
				else if ${p`j'_endog}!=1 {
					if `j'!=$nprice {
						replace ynp`j'_inst1 = y_inst*np`j'
						}
					}
			if ${p`j'_endog}==0 {
				if `j'!=$nprice {
					tempvar crenpy										
					bysort `crelev': egen `crenpy' = mean(ynp`j'_inst1)	
					quietly replace crepy`j'_inst1=`crenpy'
					}
				}
			if ${p`j'_endog}==1 {
				tempvar crenpy										
				bysort `crelev': egen `crenpy' = mean(ynp`j'_inst1)	
				quietly replace crepy`j'_inst1=`crenpy'
				tempvar crenpy
				bysort `crelev': egen `crenpy' = mean(ynp`j'_inst2)	
				quietly replace crepy`j'_inst2=`crenpy'
				tempvar crenpy
				bysort `crelev': egen `crenpy' = mean(ynp`j'_inst3)	
				quietly replace crepy`j'_inst3=`crenpy'
				}
			}

		forvalues j=1(1)$ndem {
			replace yz`j'_inst=y_inst*z`j'
			} 
		

		*with nice instrument in hand, run three stage least squares on the model, and then iterate to convergence
		replace y_old=y
		replace y_change=0
		scalar iter=0
		scalar crit_test=1
		while crit_test>$conv_crit & iter<`maxiter' {
			scalar iter=iter+1
			quietly reg3 $eqlist [aweight=wt_cur], constr($conlist) endog($endogvarlist) exog($exogvarlist)
			
			if (iter>1) {		
				matrix params_old=params
				}
			matrix params=e(b)
			quietly replace pAp=0
			quietly replace pBp=0
			quietly replace y_old=y
			quietly replace y_backup=y

			*predict with y=1
			*generate rhs vars,interactions with y=1
			forvalues j=1(1)$npowers {
				quietly replace y`j'=1
				} 
			forvalues j=1(1)$neq {
				quietly replace ynp`j'=np`j'
				tempvar crenpy 						
				bysort `crelev': egen `crenpy' = mean(np`j')
				quietly replace crepy`j'=`crenpy'
				} 
				
			forvalues j=1(1)$ndem {
				quietly replace yz`j'=z`j'
				} 

			*generate predicted values
			forvalues j=1(1)$neq {
				quietly predict s`j'hat_y1, equation(s`j')		
				}

			*set all p, pz, py to zero
			foreach yvar in $nplist $ynplist $npzlist $crepylist $creplist {  
				quietly replace `yvar'=0
				} 

			forvalues j=1(1)$neq {
				quietly predict s`j'hat_y1_p0, equation(s`j')		
				}
			
			*refresh p,pz
			forvalues j=1(1)$neq {
				quietly replace np`j'=np`j'_backup
				tempvar crenp						
				bysort `crelev': egen `crenp' = mean(np`j'_backup)
				quietly replace crep`j' = `crenp'
				forvalues k=1(1)$ndem {
					quietly replace np`j'z`k'=np`j'_backup*z`k'		 
					}
				}
			
			*generate rhs vars,interactions with y=0
			foreach yvar in $ylist $ynplist $yzlist $crepylist {  
				quietly replace `yvar'=0
				} 

				*generate predicted values
			forvalues j=1(1)$neq {
				quietly predict s`j'hat_y0, equation(s`j')		
				}

			*set all p, pz, py to zero
			foreach yvar in $nplist $ynplist $npzlist $crepylist $creplist {  
				quietly replace `yvar'=0
				} 
			forvalues j=1(1)$neq {
				quietly predict s`j'hat_y0_p0, equation(s`j')		
				}

			*refresh p only
			forvalues j=1(1)$neq {
				quietly replace np`j'=np`j'_backup
				}
			
			*fill in pAp and pBp
			forvalues j=1(1)$neq {
				quietly replace Ap`j'=s`j'hat_y0-s`j'hat_y0_p0
				quietly replace pAp=pAp+np`j'*Ap`j'	
				quietly replace Bp`j'=(s`j'hat_y1-s`j'hat_y1_p0)-(s`j'hat_y0-s`j'hat_y0_p0)
				quietly replace pBp=pBp+np`j'*Bp`j'	
				quietly drop s`j'hat_y0 s`j'hat_y0_p0 s`j'hat_y1 s`j'hat_y1_p0
				}

			*round pAp and pBp to the nearest millionth, for easier checking
			quietly replace pAp=int(1000000*pAp+0.5)/1000000
			quietly replace pBp=int(1000000*pBp+0.5)/1000000

			*recalculate y,yz,py,pz
			quietly replace y=(y_stone+0.5*pAp)/(1-0.5*pBp)
			forvalues j=1(1)$npowers {
				quietly replace y`j'=y^`j'
				}
			forvalues j=1(1)$ndem {
				quietly replace yz`j'=y*z`j'
				} 
				
			*refresh py,pz
			forvalues j=1(1)$neq {
				quietly replace ynp`j'=y*np`j'_backup
				tempvar crenpy 						
				bysort `crelev': egen `crenpy' = mean(ynp`j')
				quietly replace crepy`j'=`crenpy'
				tempvar crenp 						
				bysort `crelev': egen `crenp' = mean(np`j'_backup)
				quietly replace crep`j'=`crenp'
				forvalues k=1(1)$ndem {
					quietly replace np`j'z`k'=np`j'_backup*z`k'		 
					}
				}

			if (iter>1 & conv_param==1) {		
				matrix params_change=(params-params_old)
				matrix crit_test_mat=(params_change*(params_change'))
				svmat crit_test_mat, names(temp)
				scalar crit_test=temp
				drop temp
				}
			quietly replace y_change=abs(y-y_old)
			quietly summ y_change
			if(conv_y==1) {
				scalar crit_test=r(max)
				}
			display "iteration " iter 
			scalar list crit_test 
			summ y_change y_stone y y_old pAp pBp
			}
		
		
		*note that reported standard errors are wrong for iterated estimates
		reg3 $eqlist [aweight=wt_cur], constr($conlist) endog($endogvarlist) exog($exogvarlist)  

		mat params_`b' = e(b)
		mat params_covar_`b' = e(V)
		
		mat params_all = nullmat(params_all)\params_`b'
		}

		else if `didnotconverge'==1 {
			loc tot_didnotconverge = `tot_didnotconverge'+1
			}
			
	mat define iter_dnc_`b' = (`b',`didnotconverge')
	mat iter_dnc = nullmat(iter_dnc)\iter_dnc_`b'
	}		
		
display "Out of total iterations `bootreps', `tot_didnotconverge' did not converge."
mat list iter_dnc
mat list params_all

clear
svmat params_all, names(eqcol) 
save "${OUT}\reg3_iterated_params_boot_all_SL_food${p1_endog}good${p2_endog}service${p3_endog}_np${npowers}", replace


** Find mean param vector

loc K = colsof(params_all)
loc B = rowsof(params_all) /*only count iterations that converged*/
mat params_tot = J(1,`K',0)
forvalues b=1/`B' {
	mat params_tot = params_tot + params_all[`b',1...]
	}
mat params_mean = J(1,`K',.)
forvalues k = 1/`K' {
	mat params_mean[1,`k']=(1/`B')*params_tot[1,`k']
	}

clear
svmat params_mean

save "${OUT}\reg3_iterated_params_boot_mean_SL_food${p1_endog}good${p2_endog}service${p3_endog}_np${npowers}", replace


** Find var covar matrix
*https://cameron.econ.ucdavis.edu/research/Cameron_Miller_JHR_2015_February.pdf - look at F p11.

mat vce_tot=J(`K',`K',0)
forvalues b=1/`B' {
	mat vce_`b' = (params_all[`b',1...]-params_mean)'*(params_all[`b',1...]-params_mean)
	mat vce_tot = vce_tot + vce_`b'
	}

mat vce_boot=J(`K',`K',.)
forvalues i = 1/`K' {
	forvalues j = 1/`K' {
		mat vce_boot[`i',`j']=(1/(`B'-1))*vce_tot[`i',`j']
		}
	}

clear
svmat vce_boot

save "${OUT}\reg3_iterated_vce_boot_SL_food${p1_endog}good${p2_endog}service${p3_endog}_np${npowers}", replace


